Coulomb gap, Coulomb blockade, and dynamic activation energy in frustrated 

single-electron arrays 
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We have used modern supercomputer facilities to carry 
out extensive numerical simulations of statistical properties 
of ID and 2D arrays of single-electron islands with ran- 
dom background charges, in the limit of small island self- 
capacitance. In particular, the spectrum of single-electron 
addition energies shows a clear Coulomb gap that, in 2D ar- 
rays, obeys the Efros-Shklovskii theory modified for the spe- 
cific electron-electron interaction law. The Coulomb block- 
ade threshold voltage statistics for ID arrays is very broad, 
with r.m.s. width 8Vt growing as (Vt) oc jV 1//2 with the array 
size N. On the contrary, in square 2D arrays of large size 
the distribution around (Vt) oc N becomes relatively narrow 
(SVt/{Vt) oc l/N), and the dc I-V curves are virtually uni- 
versal. At low voltages, the slope Go(T) of I-V curves obeys 
the Arrhenius law. The corresponding activation energy Uo 
grows only slowly with iV and is considerably lower than the 
formally calculated "lowest pass" energy E ma x of the poten- 
tial profile, thus indicating the profile "softness". 
PACS numbers: 73.23.Hk, 73.40.Rw, 85.35.Gv 



I. INTRODUCTION 

Arrays of small conducting islands separated by tun- 
nel junctions are one of the simplest, and most important 
systems that exhibit the Coulomb blockade and the re- 
lated effects of correlated single-electron tunneling - see, 
e.g., Refs. 1, 2. Properties of uniform arrays, at least 
within the framework of the "orthodox" theory of single- 
electron tunneling [1,2], for both ID [3-6] and 2D [6-9] 
geometries, may be readily interpreted in terms of motion 
of single-electron solitons and anti-solitons with linear 
size (in terms of island number) M max[l, (C/Co) 1//2 ], 
where Co is the island self-capacitance, and C the capac- 
itance between the adjacent islands. 

However, in practical arrays the apparently unavoid- 
able charged impurities induce random shifts of the is- 
land electrochemical potentials, that are equivalent to 
random "background charges" Qj of the islands, spread 
uniformly within the range [— e/2, +e/2] [10] - for a dis- 
cussion see, e.g., Sec. V.C of Rcf. 11. The soliton 
language is inconvenient for the discussion of this case, 
and very little had been known about properties of such 
"completely frustrated" arrays. In the pioneering work 
[12], Middleton and Wingreen have shown that the aver- 
age Coulomb blockade threshold voltage V t of frustrated 
ID arrays with Co 3> C (i.e., M ~ 1) grows linearly with 



the array size A, while the r.m.s. width, 8Vt, of its dis- 
tribution scales as A 1 / 2 . Melsen et al. [13] found that for 
ID arrays with small island capacitance (M 3> N), (Vt) 
is proportional to A 1 / 2 . Independently, Korotkov [14] 
got the same result and also showed that in the former 
limit the relative width of the Vt distribution, 5Vt/(Vt), 
tends to a constant as N — > oo. 

Even less was known about 2D arrays. For the simplest 
case of large island capacitances Co S> C, Middleton and 
Wingreen [12] have found that at N —* oo, (Vt) oc A 
and SV t /(V t ) oc A~ 2 / 3 (ln TV) 1 / 2 . To our knowledge, all 
calculations for a finite soliton size (C ~ Co) [15,16] (as 
well as some calculations for ID arrays [17,18]) have been 
of illustrative character only. 

The goal of this work was to establish the basic proper- 
ties of the frustrated ID and 2D arrays, with an emphasis 
on the statistics of their major parameters including not 
only Vt, but also the single-electron addition energy E 
and the effective activation energy Uq that determines the 
low-temperature linear conductance of the array. For the 
sake of simplicity we have studied the most interesting 
case of small island self-capacitance (Co <C C, C/A 2 ), i.e. 
of a large soliton size (M 3> A) . In contrast to such static 
parameters as V t and E, the calculation of Uq requires 
a certain model of single-electron transport; for that we 
have used the "orthodox" theory of single-electron tun- 
neling [1], that is qualitatively valid for metallic islands 
of a not very small size (practically, above ~ lnm - see 
Fig. 2 in Ref. 11). 

In order to reach reasonable accuracy of statistical 
calculations, this work has required substantial super- 
computer resources that had become available for the 
academic solid state physics community only recently. 
Our algorithms were based on the well-known Monte 
Carlo code MOSES 1.2 that is available on the Web at 
http: / /hana. physics. sunysb.edu/set/software/index. html. 



II. DENSITY OF STATES AND COULOMB GAP 

In the absence of voltage across a frustrated array, its 
main property is the statistics of single-electron addition 
energies E of its islands. It might seem that the uniform 
distribution of Qj enables a straightforward calculation 
of the statistics, using the linearity of array electrostatics: 
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where Ek is the single-electron addition energy of fc-th 
island, i.e. the energy of a single-electron soliton cen- 
tered at that island, E^ = e 2 [C~ 1 ]hk/% is its determin- 
istic (and well known [3-7]) value in the absence of back- 
ground charges, while C _1 is the reciprocal capacitance 
matrix of the array. However, due to the random charac- 
ter of Qj , the initial distribution of these charges is gen- 
erally unstable with respect to single-electron tunneling 
between the islands. Only after a series of such single- 
electron tunneling events, the system reaches a state with 
charges = Qj + enj, with some distribution {rij} of 
additional charges, corresponding to the global minimum 
of electrostatic energy [12]. 

In experiment, such "annealing" happens automati- 
cally and is typically very fast (much faster than the typ- 
ical experiment time). However, the faithful reproduc- 
tion of this process in computer simulations requires spe- 
cial attention, because the simple relaxation within the 
framework of the orthodox theory at T = would typi- 
cally result in the system trapped in one of the metastable 
"Coulomb glass" states. Generally, the annealing may be 
carried out by increasing temperature T (that affects the 
tunneling rates [1]) and then gradually decreasing it to 
zero. We have found, however, that a faster annealing 
may be achieved by the following procedure carried out 
at T = 0: for a given island fc, the single-electron addition 
energy Ek is calculated. If the energy is negative (i.e., 
the electron addition from outside is energy- favorable), 
the electron is added (n^ — > n\~ + 1), and the operation is 
repeated. If Ek is positive, the system is allowed to relax 
due to tunneling between neighboring islands. Then the 
same procedure is carried out for the addition of a hole 
to the same island (energy E' k ). This process is repeated 
sequentially for all islands, until the addition of either 
an electron or a hole to any island is energy-unfavorable 
(all Ek,E' k > 0). By definition, this is the annealed state 
of the system. Due to the electron/hole symmetry of 
single-electron tunneling within the orthodox theory [1] , 
for the annealed state the probability distributions p(E), 
p(E') of single-electron and single-hole excitations coin- 
cide. Because of the finite array size these distributions 
depend on the island position; however, they are virtu- 
ally the same for a substantial number of the islands near 
the array center. 

Figures la,b show the density of states DoS = dp/dE 
for the central islands of, respectively, ID arrays with 
(JV — 1) islands and 2D arrays with (TV — 1) x N islands for 
several values of array size N, while the insets show the 
average values of E and r.m.s. deviations, 5E, from those 
values (i.e., energy distribution width). A suppression of 
the density of states at E — > is clearly visible. This is 
essentially the "Coulomb gap", first discussed by Efros 
and Shklovskii [19] (see also Ref. [20]). In the annealed 
state, electron transfer from island k to island j can only 
lead to an increase of the total energy of the system, i.e. 
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FIG. 1. Density of states DoS= dp/dE of (a) ID and (b) 
2D arrays as functions of the energy E of single-electron ad- 
dition to one their central islands, for several values of array 
size N. Solid line in Fig. lb shows the analytical result given 
by Eq. (1). Insets show the average single-electron addition 
energy, and the r.m.s. width 8E of its distribution, as func- 
tions of N. Here, points are numerical results; solid lines show 
E = £? (0) for arrays without frustration, while the dashed the 
best fits for frustrated arrays: (a) E m = (E) = Ne 2 /8C, 
SE = 0.13iV 1/2 e 2 /G; (b) E m = (E) = (e 2 /4nC) ln(3.9iV), 
SE = 0.019(e 2 /G)ln(40iV). 
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Ej+E' k -U jk >0. (2) 

where Ujk = e 2 [C _1 ]fcj is the energy of electron-electron 
interaction. This requirement forbids a large part of con- 
figurations {rij} and hence reduces the number of pos- 
sible states, especially at small energies, just as in hop- 
ping systems [19,20]. The only difference is quantita- 
tive: the electron-electron interaction law in the arrays is 
different from the simple unscreened Coulomb potential 
Ujk = e 2 /e|r J — r^| used for the description of hopping 
conductors. 

For example, in 2D arrays of large size N ^> 1 (but 
still N -c M) the interaction is close to 

U jk = ( e 2 /27rC)ln(0.757V/|r, - r fe |), (3) 

for 1 <C Tj,Tk <C N, i.e. for islands close enough to 
the array center. Repeating the well-known arguments 
[19,20] for this interaction law, we get the low-energy 
density of states 

DoS(E) « 13.5(C/e 2 7V 2 ) exp(25C£y e 2 ). (4) 

This dependence (shown by solid line in Fig. lb) is 
in a very reasonable agreement with the numerical re- 
sults, without any fitting parameters; in particular the 
predicted dependence DoS(O) oc iV~ 2 is satisfied sur- 
prisingly well, given the approximate character of scal- 
ing arguments and of Eq. (3). On the other hand, the 
similarity between the single-electron arrays and hopping 
conductors is not quite surprising, because the state of 
global energy minimum achieved at annealing does not 
depend on whether long electron hops are permitted (at 
hopping) or forbidden (at single-electron tunneling be- 
tween adjacent islands). 

For ID arrays, however, the similar argumentation 
does not work at low excitation energies, because in this 
case the electron-electron interaction potential [3] 

Ujk = (e 2 /2C)(N/2-\vj-v k \) (5) 

does not approach zero within the range of validity (1 -C 
fj, ft •C N) of Eq. (5), and the usual arguments ignoring 
the system boundaries [19,20] fail. Notice that in this 
case the Coulomb gap is almost "hard": for large N, 
DoS is exponentially low within a finite low-energy range 
of the order of E^ oc Ne 2 /C. 

Returning to Fig. 1, we should notice that the posi- 
tion of the DoS maximum, as well as the average single- 
electron addition energy (E) , are both very close to that 
(J<;( )) for arrays without random background charge. 
(For ID case, = Ne 2 /8C [3], while for 2D arrays, 
E(°) « (e 2 /47rC)ln(aiV), with a = 3.9 ± 0.1). This 
means, in particular, that for ID arrays even the ul- 
timately large charge frustration does not change the 
single-electron addition energies considerably. 



III. DC I-V CURVES AND COULOMB 
BLOCKADE 

Figures 2a, b show several dc I-V curves of the frus- 
trated arrays for T = 0. One can see that for 2D arrays 
of large size the statistical variations of the curves are 
really small, and in this sense the curves are virtually 
universal. (They do not depend on whether the arrays 
had been annealed before the dc voltage application, be- 
cause the motion of single-electron solitons at V > Vt 
performs effective "dynamic annealing" of the systems.) 
It means that the Coulomb blockade threshold voltage 
Vt (at which current vanishes if T = 0) also has a rel- 
atively narrow statistical distribution. This distribution 
is shown in Fig. 3b; in contrast with the case C -C Co 
[12], its r.m.s. width, expressed in units e/C, is virtually 
constant while (Vt) grows as N, just like the threshold 
(y^ )) for the case of zero background charge. (The ratio 
(V t )/V t {0) is close to 0.3.) 

On the contrary, in ID arrays the statistical distribu- 
tion of dc I-V curves (Fig. 2a), and hence the threshold 
voltages (Fig. 3a), is much broader: the V t distribution 
width grows with N almost as fast as (Vt) (as N 1 / 2 ), 
so that the relative width SVt/(Vt) decreases very slowly 
(from 0.5 to 0.4 between TV = 2 to 100). Thus we have 
confirmed (with somewhat better accuracy) the prior ID 
results [13,14]. 

As a by-product, we have found that, just like in the 
opposite case C C Co [12], the shape of dc I-V curves 
at V ^ V t is close to a power law I oc (V — Vt)^ , with 
C = 1.0 ± 0.1 in ID and ( = 1.7 ± 0.1 in 2D. 

IV. LINEAR CONDUCTANCE AND DYNAMIC 
ACTIVATION ENERGY 

Even small thermal fluctuations lift the Coulomb 
blockade [1,2] and ensure a low but finite conductance 
G = I/V of the array at low applied voltages \V\ <C V t . 
Figure 4 shows that in frustrated arrays this conductance 
is also random, with ID arrays again much more irrepro- 
ducible than square 2D arrays of the same length. Figure 
5 shows that, in contrast to systems with variable-range 
hopping [20], at low temperatures the conductance fol- 
lows the Arrhenius law 

G oc exp(—Uo/fcBT), (6) 

with the dynamic activation energy Uo depending on the 
random background charge distribution. 

Figures 6a,b show typical histograms of the statisti- 
cal distribution of Uo, while insets show the dependence 
of its first two momenta on the array size. In 2D ar- 
rays (Fig. 6b), neither of the dependences is very strong: 
(U ) « 0.09e 2 /C, SUo « 0.03e 2 /C. Notice that in 2D ar- 
rays without frustration (Qk = 0), the activation energy 
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FIG. 2. DC I-V curves of (a) ID and (b) 2D ar- 
rays, both without background charge (solid lines) and 
with completely random background charge (dash-dotted 
lines). Dashed lines show the far asymptotes of the curves 
(a) ID [3, 7], V{I) -> N[IR+sgn(I)e/2C] and (b) 2D 
V(I) — > [IR+sgn(I)Ne/4C], that are independent of the 
background charge. R is the tunnel resistance of the single 
junction between the adjacent islands. 




FIG. 3. Statistical distributions of the low-temperature 
threshold voltage Vt for (a) ID and (b) 2D arrays of various 
size. The insets show the values of Vt without background 
charge (solid curves), and the average Vt (dashed curves) and 
the r.m.s. width of their statistical distributions (dot-dashed 
curves) for frustrated arrays. In inset (a), the dashed line 
shows the phenomenological dependence (Vt) ~ 0.6(N ll/2 — 1) 
suggested in Ref. 13; the other lines are just guides for the 
eye. 
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FIG. 4. Typical dc I-V curves of several (a) ID and (b) 
2D arrays for low but finite temperatures, plotted in log-log 
scale. The curves for T = are shown with dashed lines for 
comparison.) The plots clearly show the linear conductance 
in the arrays at low voltages. 
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FIG. 5. The low-voltage conductance of ID (dashed lines) 
and 2D (solid lines) arrays as a function of inverse temper- 
ature, for several random distributions of the background 
charge. At low temperatures the plots are linear, indicative 
of the Arrhenius law: G oc exp(— Uo/ksT). 

' is also virtually independent of N: ~ 0.23e 2 /C. 

On the contrary, in ID arrays the average activation 
energy grows, albeit weakly, as N is increased (Fig. 
6a). It remains, however, much lower than its level 
(U ) « = Ne 2 /8C in arrays without frustration. 

It is also important to notice that the potential profile 
of strongly frustrated single-electron arrays is " soft" , in 
the following sense. If we plot the single-electron addition 
energy of an annealed ID array with V = as a function 
of the island number k, we get a potential profile Ek- 
It is tempting to assume that the activation energy Uq is 
just the maximum, E max , of this profile, because in order 
to carry an electron through the array at small V, ther- 
mal fluctuations should overcome this threshold. (For 
2D case, E max is determined as the lowest maximum of 
all paths connecting the opposite electrodes.) However, 
our numerical experiments have shown that this assump- 
tion is not true: Uq is always lower than E max in any 
arrays (especially those with frustration). The simple 
interpretation of this fact is that when thermal fluctua- 
tions induce the transfer of an addition electron over the 
potential profile of a frustrated array, its Coulomb field 
causes shifts of the neighboring charges and thus deforms 
the profile, reducing its maximum. 

One more important fact is that in frustrated arrays 
of any dimensionality the statistical distribution of Uo is 
very broad: 5U /(U ) oc 0.2 - 0.3 for 10 < N < 40. 



V. COMPARISON WITH EXPERIMENT 

Until recently, most experimental studies have been 
carried out on systems with random island size and loca- 
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FIG. 6. Histograms of the dynamic activation energy Uo 
for (a) ID and (b) 2D frustrated arrays. (For clarity, they 
are shown only for the largest and smallest values of N stud- 
ied in this work.) The insets show the values of Uo without 
background charge (solid curves), and the average Uo (dashed 
curves) and the r.m.s. width of their statistical distributions 
(dot-dashed curves) for frustrated arrays. The calculated val- 
ues are shown by dots; curves are only guides for the eye. 



tion - see, e.g., Refs. 21-23. In this case, the tunnel con- 
ductance between islands has exponentially broad statis- 
tical distribution, and this randomness dominates that 
created by the background charge frustration. However, 
there is a growing number of experiments [25-29] with 2D 
arrays of virtually similar and well-ordered islands, e.g. 
gold [25,27-29] or cobalt [26] clusters separated by ap- 
parently similar organic shells. (One can add the earlier 
work [24], where ID and 2D arrays have been fabricated 
by nanolithography.) For such arrays, the charge frus- 
tration should determine the property statistics, and the 
comparison with our results might be meaningful. Un- 
fortunately, for all these arrays the ratio Co/C is of the 
order of one, excluding the possibility of direct quanti- 
tative comparison of our results with the experimental 
data reported in these publications. (Our future plans 
include the extension of our calculations to this case.) 

However, several features of the experimental results 
are in a qualitative agreement with our calculations. In 
particular, 

- for large 2D arrays, the dc I-V curves are virtually 
reproducible (in particular, SVt <C (14)); 

- the general shape of these curves is rather close to 
that shown in Fig. 2b; 

- the linear conductance does obey the Arrhenius law 
(6) at low temperatures. 



VI. DISCUSSION 

The results discussed above show that the average val- 
ues and distribution widths of activation energy Uo, the 
Coulomb gap, and the Coulomb blockade threshold V t 
of charge-frustrated arrays substantially differ from each 
other, so that all these three notions should be clearly 
distinguished. (This has not always been done in earlier 
publications.) 

These results have significant practical implications. 
Two important applications for multi-junction arrays 
(sometimes called "multiple tunnel junctions", or MTJ) 
in future single-electron circuits have been suggested: 

(i) as a replacement for the double junction in single- 
electron transistors, in order to raise the single-electron 
addition energy and hence increase the operation tem- 
perature range [21,30]; and 

(ii) as a circuit element with "sub-electron" (quasi- 
continuous) charge transfer, capable of leaking to ground 
the random background charge of single-electron islands 
[11,31]. 

For application (i), the Coulomb blockade threshold of 
the array should be well reproducible (device-to-device), 
and controllable by voltage V g applied to an external 
gate electrode. Figure 3b shows that in frustrated 2D 
(but not ID!) arrays with large N the former condition 
is actually fulfilled (SVt *C (Vt)); however, the second 
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condition cannot be met. In fact, the action of exter- 
nal gate is always equivalent [1] to insertion of an addi- 
tional charge {Qk)g oc V g into each island of the array. 
If the array is completely frustrated, its random back- 
ground charges (mod e) are uniformly distributed within 
the interval [— e/2, +e/2], and the addition of additional 
charges (Qk) g does not change this statistics, and hence 
(Vt). 

The minimum requirement for application (ii) is to 
have Coulomb blockade deeply suppressed [31]: Vt <C 
e/C e , where C e is the effective capacitance of the island 
the array would shunt. Figure 3 shows, however, that for 
both ID and 2D arrays with N 3> 1 (this condition is 
necessary for reproducibility, i.e. for having SVt <C (Vt)) 
the average value of V t remains much larger than e/C. 
This means that capacitance C (and consequently the is- 
land size) of the array islands should be much larger than 
those of the island the array should serve. For practical 
integrated circuits this is hardly a good solution. (For 
discrete devices and fundamental experiments, however, 
the shunting idea may work.) 

Notice that we have explored arrays with ultimately 
low self-capacitance Co of array islands. Elementary ar- 
guments show that for both applications discussed above 
the self-capacitance Co would have the detrimental effect; 
hence our negative conclusion cannot be affected by this 
limitation. 
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